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Abstract 

The problem of hydrogen neutrals interacting with the heliospheric bow shock region has received 
considerable attention recently, motivated primarily by the hope that the Voyager spacecraft may 
soon encounter the first of the heliospheric boundaries. The complexity of the charge exchange 
interactions has limited our analytic understanding so far. In this work we develop a semi-analytic 
one-dimensional model based upon a double expansion technique to investigate the self-consistent 
interaction of shock with neutrals. The underlying method uses the Boltzmann transport equation 
that describes neutral transport, and is coupled to plasma protons predominantly through charge 
exchange processes. The gyrotropic distribution function of the neutrals is first expanded in terms 
of Legendre polynomials followed by associated Laguerre harmonics. The resulting set of equations 
can be cast into a square matrix whose eigenvalues depend upon sources and the plasma distribu- 
tion. The interstellar protons, in a high (3 limit, are described by a hydrodynamic fiuid with source 
terms that couple neutrals and plasma and further modify both distributions in a self-consistent 
manner. 
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I. INTRODUCTION 



The interaction of the supersonic solar wind with the local interstellar medium (LISM) is 
a problem of great topicality, as well as promising fundamental insights for stellar systems 
interacting with the interstellar envirnoment in general. The LISM is essentially a partially 
ionized gas, and the coupling of atoms to the solar wind and LISM plasma greatly com- 
plicates the interaction [1]. Such interactions have been been investigated using numerical 
models such as a multifluid description where plasma and neutrals are treated hydrodynam- 



icallv H, o. by t.eatmg neutral. Hnetically and the pla.na a. an unn.agnet.ed hyd.ofluW 
[3|, IJ] . The fluid description of neutrals in the outer heliosphere can however be criticized 
because their mean free paths are typically comparable to the size of the heliosphere. The 
virtually coUisionless neutrals therefore cannot equilibrate thermally on heliospheric length- 
scales which essentially prevents the use of a single fluid neutral description. 

The modeling the solar wind-LISM interaction has been based almost exclusively as 
coupled multi-dimensional models which include neutrals more-or-less self-consistently. By 
contrast, very little analytical work on the plasma- neutral coupling has accompanied the 
large-scale modeling, and this is hampering our understanding of the very complicated in- 
teraction. Here we develop a simple one- dimensional (ID) transport equation for neutral 
atoms and investigate self-consistently their influence on the structure and properties of 
hydrodynamic shocks. Since the neutrals are coupled to the plasma via charge exchange, 
we might expect the plasma to be cooled and the neutrals to be heated downstream of the 
shock. However some heated neutrals can propagate back upstream, where they can undergo 
secondary charge exchange and so heat the upstream unshocked fluid. Clearly one can ex- 
pect neutral atoms interacting with a shock wave to modify both the structure and possibly 
the character of a shock. While obviously of interest to the solar wind-LISM interaction, 
this work will have interesting implications for any shock embedded in a partially ionized 
medium, ranging from supernova shocks (SNR) to cometary bow shocks. 

Beginning with the Boltzmann transport equation for the neutrals, we employ a double 
expansion method to compute the neutral particle distribution. The production and loss of 
neutrals are computed similarly and couple to the plasma flow for which a hydrodynamical 
description is adopted. 
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II. TRANSPORT OF NEUTRALS 



The transport of neutrals across the plasma shock is governed by the time dependent 
Boltzmann transport equation describing the evolution of the neutral distribution function 
/ = /j^(x, v,t) where x, v,t denote the position, velocity and time co-ordinates; thus 

df F 

^ + V ■ V/ + - ■ Vv/ = P(x, V, t) - L(x, V, t) (1) 
ot m 

where P and L describe the production and loss of neutrals, mainly due to charge exchange, 
and are given by ^ 

P(x, v,t) = /p(x, v,t) j /(x, v',t)|v - vVexC^^v'; 

L(x,V,t) = /(x,V,t) j /(x,v',t)|v-vVexC?V. 

Here fp is the plasma distribution function, v' is the speed of neutrals and dex is the charge 
exchange (between neutrals and plasma protons) cross-section. F represents forces experi- 
enced by the neutrals due to gravity and raditation pressure. Let us now assume that the 
neutral particles are gyrotropic and stationary in a frame of the moving plasma fluid. The 
neutral speed can then be decomposed according to v = v' + u, where u = Up is bulk flow 
velocity of the plasma. Upon using spherical co-ordinates, the neutral velocity can be trans- 
formed as V = f sin 9 cos 0e^. + v sin 6 sin (pCy + v cos Ocz- We define cos 6 = fi. The transport 
equation describing the evolution of a neutral gyrotropic distribution function along the 
mean flow direction is given by 



1 — duz 
^ 2 ^ 



df ^l-fi' 



dv 2 



— + u- V m^ — V ■ u 

mv V \ot j 2 



2F, 2 d \ ^ ^ duz 

— + U- VM2 + yUV-U - 

mv V \ot / oz 



djji 



P(x,v,t)-L(x,v,t), (2) 



where Uz is the z component of the plasma velocity. The neutral transport Boltzmann 
Eq. ([2]) is rather complicated by itself and is coupled to the plasma evolution in two ways. 
Firstly, the spatial and temporal evolution of plasma velocity modifies the neutral distribu- 
tion. Secondly, the sources, i.e. the production and the loss of neutrals by charge exchange, 
depend upon the relative speed of the neutrals and the solar wind plasma proton species. 
The absolute magnitude of this relative speed can be expressed in terms of Legendre poly- 
nomials as |v — v'l = Yl'^=oOjn{'ViV')Pn{cos^), whcrc the coefficients a„ are to be determined 



by orthogonality. Using the addition theorem of spherical harmonics for the orthogonal 
Legendre polynomials and carrying out the angular integration between and 2tt allows the 
production and loss terms for the neutral species to be written as 

/■oo r + 1 °° 

P(x, V, t) = 2nfpiv) / / v'''fiv')a,^ ^ a^iv, v') PM Pnif-iW dv' (3) 

POD r + l 

L(x,v,t) = 27r/(t;) / / v'^ fp{v')a,,Y. ^n{v,v')PMPnif-i')df^'dv' (4) 

•^0 n=0 

The charge exchange parameter ciex has a logarithmically weak dependence on the relative 
speed of the neutrals and plasma. It is therefore assumed to be constant throughout our 
analysis. 



III. SOLUTION BY EXPANSION 



The multidimensionality of the neutral distribution function, described by Eq. ([2]), poses 
severe difficulties in its analytic solution. We are therefore primarily concerned with reduc- 
ing the dimension of the Eq. ([2]). For this purpose, the analytic approach developed by 
Zank et al has been used to reduce, first, the jj dependence of the neutral particle dis- 
tribution function. The distribution function is expanded in terms of Legendre polynomials 
as f{z, fi, V, t) = I]^o(2"- + ^)Pn{fJ')fn{z, V, t) . With the help of the orthogonality condition 
on the Legendre polynomials J^^ Pn{x)Pmix) = 5m,n/(2^ + 1)) where 6m,n is the Dirac delta 
function equals unity for m = n and zero otherwise, the nth component of the evolution of 
the neutral distribution function fn{z,v,t) is 

1 
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where = inp^/nh^Y/"^ is a dimensionless parameter and depends upon the ratio of equi- 
hbrium plasma {up^) and neutral (?t-/io) densities, M = Uzo/vth is plasma Mach number 
with as upstream shock speed and Vth is thermal velocity of plasma protons. In 
the non-dimensional Eq. iQ, we normalize the distribution function, length-scales, time 
and velocity of neutrals respectively as / = f /np^n'^^"^ /v^i^, z = zaf.^{np^nh^Y/'^ ^ t = 
t/vthO'cxinpf^rihoy^'^, V = v'/vth- We have neglected the force terms assuming that grav- 
ity and radiation pressure balance each other. The rhs of Eq. represents source terms 
due to constant charge exchange process and couples the plasma to the neutral dynamics. 
To work with this equation, the leading order harmonics, i.e. /o, /i, /2, ■ ■ ■, are good enough 
to describe a nearly complete distribution of the neutral particles as they contribute pre- 
dominantly to the entire distribution function. Since the neutral distribution function still 
depends upon position and velocity of the particles which may vary in time, we choose an- 
other expansion to eliminate the neutral velocity using fm{z, v, t) = J2m=o ^m{z, t)e~'" L'^iv) 
which doesn't lead to unphysically growing solution due to an exponentially damped weight- 
ing function in v. This results in a set of one- dimensional partial differential equations for 
Gm(-2, i) which can be cast into a square matrix form of x m order. The order of the 
matrix depends mainly upon how many harmonics are considered in the two expansions. 
The general form of the matrix equations, however, looks like 
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where dz is partial derivative w.r.t the z co-ordinate. The coefficients (C's, and rf's, whose 
actual expression are not shown here) of the coupling matrix could in principle be a compli- 
cated function of sources, the spatial and temporal evolution of the plasma species and so too 
are the eigenvalues. Here 4>n,'^n,^n, ■ ■ ■ are related to /i, /2, /s, ■ ■ ■ by means of the Laguerre 
polynomial expansion. The boundary condition (BC), which is crucial in this problem, can 
be set appropriately using the two expansions. For a typical problem, we inject Maxwellian 
neutrals, thermally equilibrated with respect to the plasma, through the left boundary, while 
a Dirichlet type-BC is set for the right boundary in the simulation. The entire system of 
equations can be integrated numerically in the presence of a high-/? (/3 being the ratio of 
pressure and magnetic energy of plasma particles) supersonic plasma for which we adopt a 
hydrodynamical description, as discussed in the subsequent section. 



IV. PLASMA DYNAMICS 

The plasma protons can be described by the usual compressible gasdynamical equations 
in which neutrals can influence the evolution of plasma proton in a self-consistent manner. 

| + V-(pu)=0, (7) 
^ (pu) + V ■ (puu) + Vp = Q„ (x, t ) , (8) 

i (^^^' + 7^) + ^ ■ (^^"'" + = 

While the plasma proton density is conserved, i.e. Eq. ([7]), the momentum and the energy, 
Eqs. (IE]) & (E]), can be changed by the sources Qm and Qe respectively. The sources further 
couple the plasma fluid with the neutral flow through the charge exchange processes. These 
can be described by transfer integrals as computed below j^, 

Qm(x,t) = J J (vh -Vp)\vH - Vp|(Tcx//f(x, VH,t)/p(x, Vp,t)(i^V/^(i^Vp, 

For a Maxwellian one-dimension plasma distribution, the above integrals reduce to the 
normalized form 

QU^, t) = ^ /, l^^v" - v,W + ^v't^l f{v')dv', (10) 
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f{v')dv\ 



(11) 



where all the symbols have their usual meanings. Equations ([7]) to flTTl) . describing the 
evolution of the plasma protons, along with the nth order matrix equation (for the neutrals), 
to form dynamically a self-consistent problem. 



V. CONCLUSION 



We have developed a self-consistent one-dimensional semi-analytic model for transport 
of the neutrals governed by the time dependent Boltzmann equation. The coupling of the 
neutrals to plasma protons, in a high (3 coUisionless interstellar flow, occurs predominantly 
through the neutral-proton charge exchange. The latter influences both the neutral and 
the plasma dynamics significantly. A double-expansion method for the truncated neutral 
Boltzmann equation was developed to seek a self-consistent solution of the neutral-plasma 
coupled system. The numerical simulation of the entire system of equations is underway 
and will be presented elsewhere. 
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